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i  a*  abstract  fMrnmmmrnom 

Enhancements  to  the  infrasound  software  tool  kit,  InfraMAP,  have  been  integrated  in  three  main  areas; 
propagation  modeling,  environmental  variability  modeling,  and  stochastic  localization  techniques.  All  new 
model  functionality  is  included  in  a  next-generation  release  of  the  tool  kit.  The  modeling  advances  improve 
propagation  predictions  and  understanding  of  environmental  effects,  and  they  have  been  used  to  evaluate  the 
localization  performance  and  confidence  bounds  of  operational  infrasonic  networks. 

New  propagation  modeling  features  include  a  high  altitude,  low  frequency  absorption  model,  a  synthetic 
waveform  generator  from  ray  tracing,  and  an  improved  Parabolic  Equation  (PE)  propagation  algorithm.  In 
environmental  variability,  a  range-dependent  spectral  gravity  wave  model  has  been  developed  that  generates 
wind  perturbation  fields  for  use  in  evaluating  propagation  variability.  Finally,  regional  infrasound  networks 
can  be  defined  and  used  to  compute  source  localizations  and  associated  areas  of  uncertainty,  based  upon  both 
measurement  data  and  modeling  results. 
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CONVERSION  TABLE 

Conversion  Factors  for  tf.S.  Customary  to  metric  (SI)  units  of  measurement. 


MULTIPLY - ►BY - ►TO  GET 

TO  GET  < - BY  < - - - -  DIVIDE 


angst  ran 

1.000  000  X  E  -10 

meters  (m) 

atmosphere  (normal) 

1,013  25  x  E  +2 

kilo  pascal  (kFa) 

bar 

1 . 000  000  X  E  +2 

kilo  pascal  (kPa) 

bam 

1.000  000  x  E  -28 

meter2  (m3) 

British  thermal  unit  (thermochemical) 

1.054  350  x  E  +3 

joule  (J) 

calorie  ( t±ermochemical ) 

4.184  000 

joule  (J) 

cal  ( t  hermochemical  / cm3 ) 

4 . 184  000  x  E  -2 

mega  jcule/m2  (HJ/m2) 

curie 

3 . 700  000  X  E  +1 

+giga  bacquerel  (GBq) 

degree  (angle) 

1.745  329  X  E  -2 

radian  (rad) 

degree  Fahrenheit 

tfc  *  (t°f  4-  459 . 67)  /1 . 8 

degree  kelvin  (K) 

electron  volt 

1.602  19  x  E  -19 

joule  (J) 

erg 

1.000  000  x  E  -7 

joule  (J) 

erg/ second 

1.000  000  x  E  -7 

watt  (W) 

foot 

3 . 048  000  X  E  -1 

meter  (m) 

foot  -pound™'  force 

1.355  818 

joule  (J) 

gallon  (U.S,  liquid) 

3.785  412  X  E  -3 

meter2  (nr3) 

inch 

2 . 540  000  x  E  -2 

meter  (m) 

jerk 

1 , 000  000  X  E  +9 

joule  (J) 

joule/ kilogram  (J/kg)  radiation  dose 

absorbed 

1.000  000 

Gray  [Gy) 

kilotons 

4.183 

tera joules 

kip  (1000  Ib£) 

4.448  222  X  E  +3 

newton  (N) 

kip/ inch2  (ksi) 

6. 894  757  X  E  4-3 

kilo  pascal  (kPa) 

ktap 

1.000  000  X  E  4-2 

newton  -  second/ m2  [N-s/ m2 ) 

micron 

1.000  000  X  E  -6 

meter  (m) 

mil 

2 . 540  000  X  E  -5 

meter  (m) 

mile  (international) 

1.609  344  X  E  4-3 

meter  (m) 

ounce 

2.834  952  X  E  -2 

kilogram  (kg) 

pound-forte  (Lbs  avoirdupois) 

4.448  222 

newton  (N) 

pound- force  inch 

1.129  848  X  E  -1 

newton-meter  (tf-m) 

pound  -  force/ inch 

1.751  268  X  E  4-2 

newt  on/me  ter  (N/m) 

pound- force/  foot2 

4.788  026  X  E  -2 

kilo  pascal  (kPa) 

pound-  force/  inch2  (psi ) 

6.894  757 

kilo  pascal  (kPa) 

pound-mass  (Ibn  avoirdupois) 

4 . 535  924  X  E  -1 

kilogram  (kg) 

Dound -ness -foot2  (moment  of  inertia) 

4.214  011  X  E  -2 

kilogram- meter2  (kg-m2) 

iT. 

I 

to 

1.601  846  X  E  4-1 

kilogram -me ter2  (kg/m3) 

rad  (radiation  dose  absorbed) 

1.000  000  X  E  -2 

**Gray  (Gy) 

roentgen 

2.579  760  X  E  -4 

coulomb/ kilogram  £C/kg) 

shake 

1.000  000  X  E  -8 

second  (s) 

sluq 

1.459  390  X  E  4-1 

kilogram  (kg) 

.j 

torr  (rrm  Bg,  0°  C) 

1.333  22  X  E  -1 

kilo  pascal  (kPa) 

•The  bacquerel  (Bq)  is  the  ST  unit  of  radioactivity;  1  Bq  *  1  event/s 
••The  Gray  (GY)  is  the  SI  unit  of  absorbed  radiation, 
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Section  1 
Summary 


InfraMAP,  Infrasonic  Modeling  of  Atmospheric  Propagation ,  is  a  software  tool  kit  that  has  been 
developed  for  use  by  researchers  and  analysts  active  in  the  study  of  infrasonic  propagation  and 
monitoring.  The  original  development  effort  is  summarized  in  a  previous  report  (Gibson  and 
Norris  2002b).  This  report  addresses  enhancements  that  have  been  developed  and  integrated  into 
new  versions  of  the  tool  kit. 

Enhancements  to  InfraMAP  have  been  integrated  in  three  main  areas:  propagation  modeling, 
environmental  variability  modeling,  and  stochastic  localization  techniques.  All  new  model 
functionality  is  included  in  a  next-generation  release  of  the  tool  kit.  The  modeling  advances 
improve  propagation  predictions  and  understanding  of  environmental  effects,  and  they  have  been 
used  to  evaluate  the  localization  performance  and  confidence  bounds  of  operational  infrasonic 
networks. 

Propagation  modeling  enhancements  include  a  high  altitude,  low  frequency  absorption  model,  a 
synthetic  waveform  generator  from  ray  tracing,  and  an  improved  Parabolic  Equation  (PE) 
propagation  algorithm.  In  addition  to  higher  fidelity  predictions,  the  new  functionality  enables 
more  flexibility  in  making  model-to-model  and  model-to-data  comparisons. 

In  environmental  variability,  a  new  environmental  variability  model  has  been  integrated.  Wind 
perturbations  due  to  gravity  waves  are  modeled  in  the  spectral  domain,  and  both  height  and 
range-dependent  gravity  wave  dependencies  are  modeled.  The  spectral  model  is  used  to  generate 
multiple  realizations  of  the  perturbed  environment,  which  in  turn  are  used  to  quantify  the 
statistical  variability  of  the  propagation. 

Finally,  a  network  performance  model  has  been  added,  The  model  supports  definition  of  regional 
infrasound  networks  consisting  of  three  of  more  sensor  stations,  Measurement  data  and 
modeling  results  across  these  networks  are  used  to  localize  sources  and  compute  associated  areas 
of  uncertainty.  These  capabilities  allow  for  straightforward  evaluation  of  the  effects  of 
measurement  error,  propagation  uncertainty,  and  sensor  configuration  on  source  localization 
accuracy. 


1 


Section  2 
introduction 


2.1  Background 

Monitoring  for  nuclear  explosions  requires  the  ability  to  detect,  localize,  and  discriminate 
nuclear  events.  For  the  Comprehensive  Nuclear-Test-Ban  Treaty  (CTBT).  this  monitoring  is 
being  implemented  over  a  global  scale  through  installation  of  60  infrasonic  stations  in  the 
International  Monitoring  System  (IMS).  Compared  to  other,  more  mature  monitoring 
technologies,  infrasound  has  a  greater  number  of  unresolved  issues,  most  of  which  involve  the 
effects  of  atmospheric  dynamics  on  propagation. 

The  standard  method  of  infrasonic  localization  is  based  on  taking  measured  back  azimuths  and 
arrival  times  and  projecting  them  back  to  a  source  position.  However,  due  to  the  temporal  and 
spatial  variability  of  the  atmosphere,  stable  paths  cannot  be  easily  predicted  for  standard  back- 
azimuth  and  time-o f-arrival  localization  algorithms.  It  is  generally  anticipated  that  the  ability  to 
identify  infrasound  phases,  and  include  the  appropriate  travel -time  and  bearing  refraction 
corrections  for  each  arrival  into  localization  procedures,  would  dramatically  improve  localization 
performance.  The  degree  to  which  these  capabilities  can  be  developed,  and  the  specifics  of  how 
they  should  be  used  are  major  topics  of  research  in  the  community. 

Another  infrasound  issue  of  interest  is  waveform  amplitude.  Geometrical  spreading,  diffraction, 
scattering,  and  absorption  all  have  potential  effects  on  amplitude  predictions.  Predictions  are 
necessary  for  determining  detection  thresholds  and  observability  of  individual  phases.  They  can 
also  be  used  in  source  discrimination,  in  particular  for  yield  estimates. 

Under  a  previous  DTRA  contract  (DSWA01-97-C-0160).  BBN  has  taken  a  significant  step  in 
addressing  the  needs  within  the  infrasound  monitoring  community  by  developing  the  InfraMAP 
software  tool  kit.  InfraMAP  ( Infrasonic  Modeling  of  Atmospheric  Propagation)  is  composed  of 
propagation  models  and  upper-atmospheric  characterizations,  integrated  to  allow  for  user- 
friendly  model  execution  and  data  visualization.  It  can  be  applied  to  predict  travel  times, 
bearings,  and  amplitudes  from  potential  event  locations  worldwide.  Temporal  and  spatial 
variability  of  the  atmosphere  is  addressed  by  allowing  range-dependent  temperature  and  winds  to 
be  incorporated  with  the  propagation  models. 

InfraMAP  integrates  three  types  of  research-grade  propagation  models  with  two  empirically 
based  atmospheric  models.  The  baseline  set  of  acoustic  propagation  models  consists  of  a  three- 
dimensional  ray  theory  model,  a  WKB  version  of  the  normal  mode  model,  and  a  continuous- 
wave,  two-dimensional  parabolic  equation  model.  The  baseline  empirical  atmospheric  models 
are  the  horizontal  wind  model  [Hedin  et  al,  1996]  and  the  extended  mass  spectrometer- 
incoherent  scatter  radar  temperature  model  [Picone  et  al,  1997].  Wind,  temperature,  and 
densities  are  modeled  from  the  surface  into  the  thermosphere  and  include  spatial,  diurnal,  and 
seasonal  effects. 
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2.2  Purpose 

Effective  interpretation  of  infrasound  signals  in  support  of  the  Comprehensive  Nuclear-Test- Ban 
Treaty  is  complicated  by  inherent  spatial  and  temporal  variability  in  the  atmosphere.  The 
evaluation  of  network  performance  for  inffasound  stations  relies  on  propagation  model 
predictions  for  specific  event  scenarios.  The  dynamic  nature  of  the  environment  and  its  impact 
on  the  predictions  must  be  an  integral  component  of  the  network  performance  evaluation. 

The  motivation  of  this  effort  is  to  develop  enhanced  propagation  models,  environmental 
variability  models,  and  stochastic  localization  techniques  in  support  of  an  integrated  network 
performance  model  for  infrasound.  The  anticipated  benefits  to  the  nuclear  explosion  monitoring 
community  are: 

Improved  modeling  capabilities  -  Propagation  modeling  enhancement  will  improve  the  fidelity 
of  the  models,  leading  to  more  accurate  estimates  of  location  and  yield.  The  model  enhancements 
will  also  facilitate  model-to-model  and  model-to-measurement  comparisons,  which  will  help 
define  each  model's  performance  and  applicability. 

Evaluation  of  IMS  network  -  Network  performance  models  will  enable  evaluation  of  the 
effectiveness  of  the  IMS  network  for  specified  scenarios  of  interest. 

Improved  infrasound  localization  -  Integration  of  environmental  variability  models  into  event- 
localization  algorithms  will  improve  event  location  and  confidence  bound  estimates.  The 
confidence  bounds  will  be  useful  in  determining  the  relative  importance  with  which  to  treat  the 
predictions. 

Support  of  synergy-based  localization  -  Confidence  bounds  attached  to  event-localization 
predictions  will  be  useful  as  weighting  factors  in  synergy-based  localization  algorithms  that 
process  predictions  from  multiple  monitoring  technologies. 


2.3  Scope 

This  effort  has  been  limited  to  integrating  improved  propagation  models,  variability  models,  and 
localization  capabilities  into  new  versions  of  the  InfraMAP  software  tool  kit.  Specifically,  the 
scope  is  limited  to: 

•  Enhancing  propagation  modeling  capabilities. 

•  Developing  spectral  environmental  variability  models. 

•  Developing  network  performance  models  that  incorporate  environmental  variability. 

•  Performing  model-to-model  and  model-to-measurement  comparisons  to  validate  models 
and  define  confidence  levels. 

•  Integrating  all  model  enhancements  and  increased  functionality  into  a  new  release  of 
the  InfraMAP  software  tool  kit. 
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2.4  Outline 


This  section  provides  an  introduction  to  the  research  effort.  The  remainder  of  the  report 
summarizes  the  technical  work  that  has  been  performed.  Section  3  addresses  methods, 
assumptions,  and  procedures  in  the  model  development  process.  Section  4  discusses  the  results 
of  the  effort  by  describing  the  software  functionality.  Studies  are  also  summarized  that  have 
been  performed  using  the  software.  The  report  ends  with  conclusions  and  recommendations. 
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Section  3 

Methods,  Assumptions  and  Procedures 


This  section  describes  the  approach  to  enhancing  the  InfraMAP  capabilities.  The  propagation 
modeling  enhancements  are  addressed  first,  followed  by  discussions  of  environmental  variability 
and  network  performance.  The  section  concludes  writh  a  summary  of  the  software  development 
methods. 


3.1  Propagation  Modeling 

Propagation  model  enhancements  were  selected  in  order  to  support  the  project  goals.  They  fall 
into  two  categories:  improvements  in  sophistication  and  efficiency,  and  improvements  to  support 
model-to-modei  and  model-to-measurement  comparisons.  Implementation  of  items  in  the  first 
category  will  naturally  lead  to  better  predictions.  That  for  the  second  category  will  facilitate 
definition  of  model  confidence  levels  and  aid  in  evaluation  of  model  applicability. 

Table  1  provides  a  summary  of  the  strengths  and  weaknesses  of  the  propagation  models  currently 
integrated  into  InfraMAP.  This  table  was  used  to  define  the  enhancements  for  integration  into 
the  three  models,  as  detailed  below. 


Table  1.  Strengths  and  weaknesses  of  baseline  propagation  models. 


Propagation  Model  Strengths  Weaknesses 

Ray 

•  Prediction/Visualization  of 
acoustic  wave  front  trajectories 

•  Three-dimensional  domain 

■  Accounts  for  moving  medium 

•  Accounts  for  range  and  azimuthal- 
dependent  environment 

•  Accounts  for  evolution  of 
environment 

*  Poor  amplitude  predictions 

*  Ray  path  approximation  breaks 
down  at  low  frequencies 

*  High  numerical  load  for  many 
ray  traces 

*  No  signal  waveform 
predictions 

Normal  Mode 

•  Prediction  of  signal  waveform 

*  Accounts  for  source  signal 
characteristics  (in  time  domain) 

•  Two-dimensional  domain 
(vertical  slice) 

•  Accounts  for  only  range- 
independent  environment 

•  Potential  for  false  reflections  at 
upper  boundary 

Parabolic  Equation  (PE) 

•  Prediction/Visualization  of 
amplitude  field  down  range  of 
source 

*  Accounts  for  range-dependent 
environment 

•  Accounts  for  source  frequency 

*  Accounts  for  acoustic  ground 
impedance 

*  Two-dimensional  domain 
(vertical  slice) 

*  Sensitivity  to  source-field 
characterization 

*  Limited  to  continuous  wave 
(CW)  sources 

1  1 
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3.1.1.  Ray  Model 


The  ray  model  integrated  into  InfraMAP  is  HARPA  (//Amiltonian  Ray  Propagation  in  the 
Atmosphere)  [Jones  et  al.,  1986],  it  is  defined  in  a  three-dimensional  domain  by  spherical 
coordinates.  HARPA  accounts  for  the  moving  medium  and  predicts  ray  path  attenuation  due  to 
viscous  and  thermal  losses.  Interaction  with  the  surface  is  handled  with  specular  reflections. 

Under  BBN's  efforts  in  CTBT  hydroacoustic  modeling,  it  was  found  that  the  difference  in  ray 
paths  traced  over  an  elliptical  earth  model  versus  a  spherical  earth  model  was  non-trivial  [Farrell 
and  LePage,  1996].  The  conclusion  was  that  an  elliptical  earth  model  should  be  used  for  ray 
traces  and  ray- based  travel -time  grids,  and  the  implementation  of  HARPA  in  InfraMAP  has  been 
enhanced  by  incorporating  an  elliptical  earth  model,  This  enhancement  enables  evaluation  of  ray 
sensitivity  to  the  earth  characterization. 

Absorption  can  have  a  significant  influence  on  detection  of  infrasound  phases.  For  example, 
rays  that  sample  the  atmosphere  above  120  km  experience  high  levels  of  absorption  and  are 
generally  not  observed  [Gibson  et  al. ,  1998].  The  absorption  model  used  with  HARPA  in 
InfraMAP  has  been  improved  by  incorporating  the  recent  mode!  developed  by  Sutherland  and 
Bass  [1996],  This  absorption  model  is  concluded  to  best  account  for  the  relevant  physics  of 
infrasound,  as  it  was  specifically  developed  with  low-frequency  propagation  in  mind. 

To  facilitate  comparison  of  ray  predictions  to  normal  mode  predictions,  a  signal  phase  model  has 
been  incorporated  with  HARPA  in  InfraMAP.  The  signal  phase  is  predicted  from  a  user-defined 
source  function  and  a  set  of  rays.  Although  absolute  predictions  of  ray  levels  are  of  limited 
application,  this  technique  has  been  used  successfully  in  the  ocean  to  associate  ray  arrivals  with 
measured  waveforms  [Norris  et  al.  1996],  Each  ray  is  used  to  superimpose,  at  the  ray  travel 
time,  the  source  waveform  weighted  by  the  ray  attenuation.  Predicted  phases  of  this  type  can  be 
directly  compared  to  normal  mode  predictions.  In  addition,  multiple  source  functions  can  be 
used  in  comparisons  with  measurements  to  help  quantify  the  source  characteristics. 

To  further  facilitate  model-to-model  comparisons,  two  environmental  options  have  been 
incorporated  in  the  ray  model.  The  first  option  reduces  the  range-dependent  environment  to  a 
range-independent  one,  i.e.,  a  single  environmental  profile,  specified  at  the  source,  receiver,  or 
averaged  along  the  source-receiver  path.  The  second  option  reduces  the  time-dependent 
environment  to  a  time-independent  one,  i.e.,  an  environmental  field  frozen  at  a  specified  time. 
These  comparison-based  options  are  desired  because  range -dependent  and  time-dependent 
environments  are  not  accounted  for  in  the  other  propagation  models  in  InfraMAP. 

3.1.2.  Parabolic  Equation  (PE)  Model 

Predicting  the  loss  of  an  infrasonic  signal  as  it  propagates  from  source  to  receiver  is  important  in 
the  study  of  network  performance.  It  helps  determine  the  receiver  station  coverage  and  detection 
threshold  for  a  given  event  scenario.  It  also  aids  in  determining  the  observability  of  acoustic 
paths  that  sample  the  thermosphere,  a  region  of  high  acoustic  absorption.  Thermospheric  paths 
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are  of  particular  interest  because,  when  analyzed  along  with  lower  atmospheric  paths,  they  can 
be  used  to  estimate  source  position. 

Loss  predictions  are  made  with  the  Parabolic  Equation  (PE)  model.  The  PE  model  originally 
integrated  into  InffaMAP  is  a  continuous- wave,  finite-difference  implementation  with  no 
absorption  model  [West  et  al.,  1992],  and  it  can  be  applied  to  a  vertical  range-dependent  slice  of 
the  atmosphere. 

The  PE  model  in  InfraMAP  has  been  improved  by  incorporating  a  more  efficient,  higher  fidelity 
split-step  algorithm  [Jensen  et  al.,  1994],  In  addition,  an  improved  source  field  realization  model 
has  been  integrated  that  more  accurately  models  an  omnidirectional  point  source.  Both  of  these 
PE  enhancements  provide  more  accurate  amplitude  predictions. 

The  PE  model  has  also  been  enhanced  by  adding  an  absorption  model.  In  the  thermospheric 
region  above  100-130  km,  absorption  increases  significantly.  The  main  effect  of  the  absorption 
is  a  reduction  in  acoustic  energy  from  this  region.  The  same  absorption  model  used  for  the  ray 
enhancements  [Sutherland  and  Bass,  1996]  has  been  integrated  here,  allowing  for  compatibility 
in  comparisons. 

Theoretical  developments  that  enable  travel  time  predictions  of  impulse  sources  from  PE  models 
have  been  evaluated  in  terms  of  feasibility  and  computational  loading.  Both  the  Time-Domain 
PE  (TDPE)  [Jensen  et  al. ,  1994]  and  the  Nonlinear  Progressive  wave  Equation  (NPE) 

[McDonald  and  Kuperman,  1987]  have  been  identified  as  promising  models  for  future 
development  and  application  in  the  study  of  infrasonic  propagation. 

3.1.3.  Normal  Mode  Model 

The  normal  mode  model  in  InfraMAP  is  a  WKB  approach  [Dighe  et  al. ,  1998;  Hunter  and 
Whitaker,  1997]  of  the  normal  mode  model  [Pierce  and  Kinney,  1976;  Pierce  et  al.,  1973;  Pierce 
and  Posey,  1970].  It  can  be  applied  to  a  range-independent  vertical  slice  of  the  atmosphere. 

False  reflections  from  the  upper  boundary  of  the  normal  mode  model  can  appear  in  the 
waveform  predictions,  interfering  with  signal  phase  interpretation.  An  evaluation  of  the  upper 
boundary  has  been  made  and  optimal  environmental  boundary  conditions  have  been  identified 
and  integrated  in  order  to  minimize  false  reflections. 


3.2  Environmental  Variability 

Environmental  properties  vary  in  both  space  and  time.  Coherent  spatial  variability  is  observed  at 
length  scales  from  meters  to  thousands  of  kilometers,  and  temporal  variability  occurs  over 
diurnal  and  seasonal  time  scales.  The  variability  in  wind  and  temperature  makes  modeling 
infrasound  propagation  difficult.  Since  the  environment  is  dynamic  and  cannot  be  measured 
over  the  vast  regions  through  which  the  infrasonic  signals  propagate,  stochastic  modeling 
methods  are  necessary  to  account  for  the  environment’s  influence  on  propagation. 
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BBN  has  taken  an  initial  step  in  this  area  by  incorporating  wind  perturbation  profiles  in  the 
InfraMAP  ray  calculations.  The  environmental  perturbation  fields  are  defined  using  a  one¬ 
dimensional  vertical  w'ave  number  spectrum  of  the  horizontal  wind,  and  individual  realizations 
of  spatial  fields  are  generated  through  a  random-phase  technique  [Peitgen  and  Saupe,  1998].  A 
Monte  Carlo  simulation  is  executed  where  multiple  rays  are  traced  through  the  sum  of  mean  and 
perturbed  profiles.  This  simulation  is  used  to  estimate  moments  of  travel-time  and  azimuthal 
deviation. 

These  modeling  capabilities  are  a  first  step  in  addressing  the  variability  issue.  However, 
improving  the  stochastic  environmental  model  is  needed  to  more  realistically  represent  the 
variability.  The  improvements  are  made  by  integrating  a  two-dimensional  horizontal  wave 
number  spectrum  model  [Gardner.  1995.  1993].  This  model  captures  the  environmental 
variability  associated  with  atmospheric  gravity  waves.  Gravity  waves  are  modeled  because  their 
spatial  scales  are  of  the  same  order  as  infrasonie  wavelengths  and  thus  they  are  suspected  of 
having  the  greatest  influence  on  propagation  variability. 

Ultimately,  by  better  defining  the  spatial  dependence  and  amplitudes  of  the  environmental 
perturbations,  better  estimates  of  the  stochastic  properties  of  the  propagation  will  be  obtained. 
These  estimates  are  an  integral  part  of  the  network  performance  predictions,  discussed  in  the 
next  section. 


3.3  Network  Performance. 

The  CTBT  monitoring  network  is  designed  to  provide  event-location  estimates  from  observed 
infrasound  signals  at  a  set  of  monitoring  stations.  The  sophistication  of  the  propagation  and 
environmental  models,  as  well  as  the  choice  of  localization  algorithm,  determines  the  accuracy 
of  these  estimates.  A  network  performance  model  has  been  developed  that  incorporates 
enhanced  propagation  modeling  and  environmental  variability.  The  performance  model  is  based 
on  predicted  event  locations  and  areas  of  uncertainty  (AOU)  for  specified  scenarios.  The  model 
includes  comprehensive  integration  of  the  uncertainties  across  a  network,  which  not  only 
provides  better  performance  predictions,  but  also  helps  quantify  the  influence  of  the  propagation 
model  enhancements  and  environmental  variability. 

The  model  uses  a  localization  algorithm  based  on  the  Maximum  Likelihood  Estimator  (MLE) 
approach  (Kay,  1993).  Other  localization  algorithms  exist,  but  it  has  been  determined  that  the 
performance  difference  between  various  algorithms  is  minimal.  Back  azimuth  and  arrival  time 
measurements  and  associated  uncertainties  are  used  in  the  localizations.  The  uncertainties  can 
include  both  measurement  error  and  uncertainty  associated  with  the  propagation.  Propagation 
variability  due  to  the  environment  is  also  incorporated  into  the  calculations,  and  estimates  of 
event  location  and  area  of  uncertainty  (AOU)  are  accomplished. 
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3.4  Software  Implementation. 

The  existing  InfraMAP  software  tool  kit  serves  as  a  development  platform  for  model 
enhancement,  creation,  and  evaluation.  New  functionality  has  been  added  to  InfraMAP  as  the 
model  development  progresses,  and  incremental  releases  of  InfraMAP  software  have  been 
distributed  to  the  user  community. 

In  the  development  of  the  software,  technical  interaction  with  end  users  has  provided  valuable 
feedback  with  respect  to  refining  the  original  design  concepts.  As  an  example,  one  user 
suggested  the  incorporation  of  the  azimuthal  equidistant  projection  to  view  topography  and 
propagation  predictions,  as  this  projection  is  useful  to  quickly  visualize  range  contours.  The 
azimuthal  equidistant  projection  is  now  an  integrated  InfraMAP  option. 

User  support  of  the  InfraMAP  software  tool  kit  has  been  provided  in  parallel  with  the  other  main 
tasks.  This  user  support  has  served  two  purposes.  First,  full  benefit  of  the  tool  kit  has  been 
ensured  by  offering  technical  help  and  customized  routines  for  specific  user  needs.  Second, 
relevant  issues  within  the  infrasonic  community  have  been  assimilated  through  technical 
interchange  and  user  feedback.  The  feedback  has  been  reviewed  for  unanticipated  user  needs, 
and.  where  appropriate,  the  focus  of  software  development  has  been  refined  to  provide  maximum 
benefit  and  applicability  to  the  infrasonic  researcher  and  analyst. 
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Section  4 

Results  and  Discussion 


This  section  details  the  results  and  accomplishments  that  have  been  made  under  this  contract. 
First,  enhancements  and  new  functionality  to  the  InfraMAP  software  tool  kit  are  reviewed.  This 
area  is  presented  with  an  initial  overview,  followed  by  detailed  discussions  in  the  areas  of 
propagation  modeling,  environmental  variability,  and  network  performance.  Second,  references 
are  given  to  the  technical  reports  and  presentations.  The  section  concludes  with  a  summary  of 
the  sensitivity  and  validation  studies. 


4,1  Overview  of  New  InfraMAP  Functionality. 

The  baseline  InfraMAP  software  tool  kit  was  released  in  2000  [Gibson  and  Norris,  2002b]. 

Many  new  major  and  minor  features  have  been  added  to  the  baseline  release.  They  generally 
involve  improving  the  graphical  user  interfaces  (GUIs)  or  improving  the  model  functionality. 
Table  2  provides  a  summary  of  new  features  added  under  this  contract.  Major  additions  to 
InfraMAP  are  discussed  in  this  section.  Full  details  concerning  InfraMAP  operation  and  options 
can  be  found  in  the  latest  User's  Guide  [Norris  and  Gibson,  2003], 


Table  2.  Overview  of  new  InfraMAP  functionality. 


General  Improvements 


Improved  Model  executables  -  All  model 
executables  are  compiled  statically,  which 
simplifies  installation  and  improves 
performance _ 

Improved  Print/Save  form  -  Data  from  any 
figure  can  now  be  saved  to  MATLAB  or  text 
files _ 

Increased  limit  on  number  of  Ray  Trace  step 
iterations  per  bounce  -  Prevents  premature  exit 
of  ray  tracing  that  can  occur  during 
perturbation  runs _ 

New  file  interface  supporting  subdirectories 
for  propagation  runfiles _ 

Run  Status  Viewer  -  During  propagation  runs, 
the  user  can  view  the  run  status  using  the  "Run 
Status”  button  on  the  Run/View  form _ 

Run  Log  Viewer  -  Run  messages  are  now 
printed  to  a  log  file  instead  of  being  displayed 
at  the  command  line.  The  user  can  view  the 
run  log  using  the  “View  Run  Log”  button  of 
the  Run/View  form 
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Updated  Ellipsoidal  Earth  Model  -  WGS84 
standard.  Spherical  earth  model  radius  set  to 
average  of  ellipsoidal  major  and  minor  axis 

Extended  Reverse  Propagation  Functionality  - 
Option  of  specifying  back  azimuth  when 
reverse  propagation  is  selected 

New  Coordinate  Display  Functionality  - 
Option  of  specifying  lat/lon  coordinates  in 
either  decimal  degrees  or  deg:min 

Main  Window  Functionality 

Improvements 

New  Map  Overlay  Functionality  -  Option  of 
overlaying  map  features  on  any  lat/lon  figure. 
Option  of  overlaying  propagation  scenario 
great  circle  path  on  any  lat/lon  figure 

Easier  Cross-Application  settings  transfer  - 
New  option  for  transferring  Propagation 

Scenario  form  data  to  Propagation  Model  form 

New  variable  options  for  MSIS  -  Molecular 
weight  and  specific  heat  ratio  computed  from 
MSIS  can  be  displayed  in  the  View 

Environment  form 

Rav  Model/Interface  Improvements 

Improved  Eigenray  finder  -  Two  passes  over 
initial  search  over  fan  of  rays  now  includes 
correction  for  azimuthal  deviation,  reducing 
the  probability  of  missing  candidate  eigenrays 

New  Ray  Synthesis  Functionality  -  Ray  path 
predictions  can  now  be  coupled  with  a  source 
waveform  to  predict  a  receiver  envelope  that 
includes  propagation  losses  due  to  spreading 
and  atmospheric  absorption.  The  user  can 
select  from  a  Pierce  blast  wave  pulse,  a 
variable  width  Gaussian  pulse  or  a  user- 
defined  source  waveform 

Sutherland  and  Bass  Absorption  Model  - 
Sutherland  and  Bass's  atmospheric  absorption 
model  can  now  be  included  in  ray  synthesis 
waveform  predictions 

Range  Independent  HW  M/MS  IS  Option 

Ellipsoidal  Earth  Option  for  the  Ray  Model 

Normal  Mode  Model/Interface 

Improvements 

Independent  Specification  of  number  of  points 
used  in  upper/lower  duct  dispersion  curve 
solutions 

Specification  of  time  window  for  averaging  of 
normal  mode  model  waveforms 

New  upper  boundary  condition  for  normal 
mode  code  that  better  models  receiver 
amplitudes  for  thermospheric  arrivals 

]  1 


Auto  scaling  option  for  amplitude  axis  of 
normal  mode  waveform  plots 

PE  Model/Interface  Improvements 

Improved  upper  boundary  condition 
(absorption  layer)  that  eliminates  false 
reflections 

Range-independent  environmental  options 
similar  to  that  for  ray  and  normal  mode 
models 

Sutherland  and  Bass  Absorption  Model  - 
Sutherland  and  Bass’s  atmospheric  absorption 
model  can  now  be  included  in  Parabolic 

Equation  propagation  predictions 

ON/OFF  button  for  Time  Evolution  of 
Environment  in  PE  Model 

New  Split-Step  Algorithm  and  Source 
Specification 

New  PE  Model  Display  Option  -  User 
specification  of  amplitude  reference  distance 

Separate  steps/wavelength  inputs  for  height 
and  range 

Propaaation  Variability  Improvements 

New  Gravity  wave  spectral  model  for  wind 
perturbation 

New  range-dependent  specification  of 
perturbation  field 

New  plots  for  Monte  Carlo  simulation  that 
show  path  of  reference  ray  and  perturbed  rays 
in  both  vertical  and  horizontal  plane 

New  option  to  fix  range  of  perturbed  ray 
tracing  -  This  option  is  useful  if  perturbed  rays 
deviate  excessively  from  reference  ray. 

Network  Performance  Improvements 

New  Source  Localization  form  that  supports 
calculation  of  source  location  and  area  of 
uncertainty  based  on  input  azimuth  and  travel 
time  data  across  selected  network 

Improved  Station  Location  and  Station 

Network  forms 

Updates  to  infrasound  station  database  to 
include  more  recent  locations  and  names 

Comprehensive  updates  to  IMS  and  prototype 
infrasound  station  locations 

4.2  Propagation  Modeling  Enhancements. 

The  three  propagation  models  (ray  trace,  normal  mode,  PE)  have  a  variety  of  capabilities  and 
environmental  interface  options.  They  are  summarized  in  Table  3.  The  PE  and  ray  trace  models 
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are  range-dependent  models,  while  the  normal  mode  model  only  propagates  through  a  range- 
independent  environment.  Three-dimensional  effects,  such  as  azimuthal  deviation,  are  only 
predicted  using  ray  tracing.  The  normal  mode  model  produces  waveform  predictions,  and  the 
ray  trace  model  can  be  used  to  generate  synthetic  waveforms  based  upon  an  assumed  source 
function.  Details  concerning  the  propagation  modeling  are  given  in  the  following  subsections 


Table  3.  Summary  of  propagation  capabilities  and  environmental  interface  options. 


Predictive 

Capabilities 

Travel  Time 

Azimuth 

Deviation 

Waveform 

Absorption 

Environemnt 

Options 

User-defined 

Profiles 

HWM/MSIS 

Range  Dep. 

HWM/MSIS 

Range  Indep. 

HWM/MSIS 

Freeze  Time 

Earth 

Curvature 

Ray  Tracing  (3-D) 

-  HARPA  (NOAA) 

-  Geometrical  acoustics 

YES 

YES 

synthetic 

YES 

YES 

YES 

YES 

YES 

YES 

Normal  Modes  (WKB) 

*  Modified  Pierce  code 

-  Multi-mode  summation 

YES 

YES 

YES 

YES 

Parabolic  Equation  (2-D) 

-  Split-step  PE 

-  Range-marching 

YES 

YES 

YES 

YES 

YES 

4.2.1  Ray  Tracing  Model. 

The  3-D  ray  tracing  code  HARPA  (Hamiltonian  Ray  Tracing  Program  for  Acoustic  Waves  in  the 
Atmosphere)  was  developed  at  NOAA  and  provides  ray  predictions  through  an  inhomogeneous 
three-dimensional  representation  of  the  atmosphere  [Jones  et  ai,  19861.  The  model  accounts  for 
vertical  and  horizontal  refraction  as  well  as  horizontal  translation  of  the  ray  path  due  to  the 
moving  medium.  Horizontal  refraction  and  translation  effects  are  important  in  that  they  can 
cause  significant  azimuth  biases  in  the  propagation  path.  HARPA  output  includes  azimuth  and 
elevation  angle,  travel  time,  and  azimuthal  deviation  at  each  point  along  a  ray  path.  The  ray- 
tracing  technique  is  capable  of  stepping  either  forward  from  a  source  or  backward  from  a 
receiver.  In  InfraMAP,  HARPA  has  been  integrated  with  an  eigenray  finder,  so  that  rays  can  be 
identified  that  connect  a  source  and  a  receiver  to  within  a  given  miss  distance. 

The  ray  tracing  model  has  been  enhanced  with  a  new  version  of  the  eigenray  finder.  The 
eigenray  algorithm  starts  by  shooting  a  fan  of  rays  over  elevation  angle  along  an  azimuth 
directed  toward  the  receiver.  A  second  fan  of  rays  is  computed  in  the  new  version.  The  fan  is 
identical  except  that  the  azimuth  for  each  ray  is  corrected  for  azimuthal  deviation.  The 
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correction  improves  the  robustness  of  the  algorithm  in  scenarios  with  long  ranges  or  strong 
crosswinds.  The  miss  distance  to  the  receiver  is  computed  for  the  rays  in  the  second  fan,  and 
rays  that  correspond  to  local  minima  in  the  miss  distance  are  identified.  These  rays  are  then  used 
as  candidate  eigenrays,  Each  candidate  ray  is  perturbed  in  both  azimuth  and  elevation  angle  until 
a  minimum  miss  distance  to  the  receiver  is  found.  The  ray  is  identified  as  an  eigenray  if  the  miss 
distance  falls  within  the  user  specified  limits. 

A  synthetic  waveform  option  has  also  been  added  for  use  with  the  ray  trace  model.  The  source 
waveform  is  first  selected,  which  can  be  defined  as  blast  wave,  Gaussian  waveform,  or  user- 
defined  waveform.  The  blast  wave  is  computed  from  the  formulation  used  in  InfraMAP’s 
normal  mode  code  (Pierce  and  Posey,  1970).  Figure  1  shows  a  blast  wrave  for  a  1  kT  source  on 
the  ground. 


time  (s) 

Figure  1.  Blast  waveform  model  for  a  1  kT  source  on  the  ground  [Pierce  and  Posey,  1970]. 


Once  the  source  waveform  is  selected,  the  synthetic  waveform  is  computed  by  convolving  the 
source  waveform  with  weighted  impulse  functions  at  the  arrival  times  of  the  selected  rays.  The 
impulse  functions  are  weighted  by  the  attenuation  loss  calculated  over  each  ray.  Figure  2 
provides  a  graphical  representation  of  the  waveform  synthesis  process. 
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Figure  2.  Graphical  flow  diagram  depicting  synthetic  waveform  process. 

4.2.2  Normal  Mode  Model. 

A  version  of  Pierce's  normal  mode  code,  modified  for  prediction  of  lower-yield  nuclear  events, 
has  been  developed  by  Los  Alamos  National  Laboratory  [Dighe  et  al,,  1998;  Hunter  and 
Whitaker,  1997],  The  model  uses  the  Wentzel-Kramer-Brillouin  (WK.B)  method  to  calculate  a 
dispersion  curve,  which  in  turn  is  used  to  calculate  a  received  waveform  at  a  given  range  and 
azimuth.  The  waveform  is  computed  from  the  summed  set  of  modal  solutions.  Only  range- 
independent  environments  are  supported,  which  can  be  based  on  source  location,  receiver 
location,  or  from  averaged  values  along  the  entire  propagation  path.  The  special  case  of  dual¬ 
duct  propagation  is  identified  and  automatically  processed,  in  which  case  the  full  waveform 
solution  is  the  sum  of  the  upper  and  lower  duct  solutions. 

The  nomial  mode  model  uses  an  isothermal  half-space  as  an  upper  boundary  condition.  In 
working  with  Los  Alamos  National  Laboratory  (Rod  Whitaker),  amplitude  predictions  for 
thermospheric  arrivals  were  observed  to  be  too  strong  when  compared  to  observations,  and 
energy  returned  from  the  upper  boundary  conditions  was  identified  to  contribute  to  this  problem. 
As  a  result,  the  isothermal  half-space  has  been  modified  to  reduce  the  strength  of  thermospheric 
arrivals. 
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4.2.3  Parabolic  Equation  Model. 

The  Parabolic  Equation  (PE)  algorithm  steps  forward  from  a  source,  characterized  by  a  starter 
field,  and  calculates  an  amplitude  field  in  height  and  range  at  a  single  frequency.  The 
atmosphere  is  defined  at  each  range  step  and  can  be  range-dependent.  The  baseline  PE  version 
used  a  finite-difference  approach.  The  model  has  been  upgraded  by  integrating  a  split-step 
Fourier  algorithm  approach  [Jensen  et  al.,  1994].  The  new  version  is  computationally  more 
efficient  and  more  accurately  models  the  source  function.  The  upper  boundary  condition  has 
been  modified  to  eliminate  false  reflections  that  could  appear  in  the  baseline  version. 

A  new  environmental  interface  option  has  also  been  added  to  the  PE  model.  PE  predictions  can 
be  made  through  an  environment  that  is  evolving  in  time.  A  mean  signal  velocity  must  be 
defined  for  a  given  scenario,  and  this  is  used  to  update  the  "local"  time  at  each  PE  marching 
solution  range,  This  time  is  then  used  in  the  retrieval  of  the  HWM/MSIS  environmental  data, 
thus  incorporating  the  environmental  time  evolution  in  the  propagation  predictions. 

4.2.4  Absorption 

An  acoustic  absorption  model  has  been  integrated  for  use  with  the  PE  and  ray  tracing  models  to 
improve  the  accuracy  of  the  amplitude  predictions.  The  absorption  model  chosen  was  developed 
specifically  for  use  with  low  frequencies  and  altitudes  up  to  160  km  [Sutherland  and  Bass,  1996]. 
It  includes  the  contributions  from  both  classical  (translation  and  diffusion)  and  relaxation 
(rotation  and  vibration)  losses.  The  model  computes  absorption  using  vertical  profiles  of 
temperature,  pressure,  and  the  concentration  of  gases  that  make  up  the  air.  All  of  these  input 
variables  are  obtained  from  the  MSISE-90  model  (Picone  et  ah,  1 997).  Figure  3  gives  an 
example  of  the  MSIS  gas  densities  that  are  used  for  computing  absorption.  See  Section  4.6.2  for 
example  absorption  coefficient  profiles. 
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density  (particies/cc) 


Figure  3,  Example  profiles  of  MSISE-90  gas  densities  used  in  the  absorption  model. 


4.3  Environmental  Variability  Enhancements 

The  baseline  environmental  variability  model  in  InfraMAP  uses  a  power-law  wind  perturbation 
spectrum,  and  it  provides  realizations  of  wind  perturbation  profiles.  This  spectrum  is  applicable 
for  small-scale  turbulence,  even  though  atmospheric  turbulence  covers  a  wade  spectrum  of 
spatial  scales.  The  dominant  source  of  variability  affecting  infrasonic  propagation  is  believed  to 
result  from  gravity  waves.  Gravity  w'aves  result  from  oscillations  of  air  parcels  displaced  by 
buoyancy  and  restored  by  gravity.  The  oscillations  have  time  scales  ranging  from  minutes  to 
tens  of  hours.  Vertical  length  scales  of  gravity  waves  are  in  the  range  of  0.1  to  10  km.  and 
horizontal  scales  can  span  from  100  to  10,000  km. 

A  significant  body  of  research  has  been  carried  out  to  define  the  spectral  character  of  gravity 
waves.  The  spectral  model  of  Gardner  [Gardner  1995,  1993]  has  been  selected  for  integration 
into  the  environmental  variability  module.  This  model  is  based  on  scale- independent  diffusive 
filtering  theory.  A  source  spectrum  is  defined  near  the  ground.  As  the  spectrum  is  propagated  up 
in  height,  attenuation  is  modeled  by  introduction  of  diffusive  damping.  The  key  spectral 
properties  are; 

*  increase  in  energy  with  height 

*  shift  towards  larger  length  scales  with  height 

*  attenuation  of  smaller  length  scales  with  height 
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The  Gardner  Spectral  model  is  evaluated  at  four  discrete  heights,  as  shown  in  Figure  4.  These 
heights  capture  the  dominant  gravity  wave  variability  from  the  troposphere  up  to  the  lower 
thermosphere.  Gravity  waves  are  not  fully  developed  below  the  troposphere,  In  the 
thermosphere,  diffusion  increases  dramatically  and  gravity  waves  are  damped  out.  In  the  figure, 
the  predicted  standard  deviations  of  the  wind  perturbations  are  listed  alongside  the  spectral 
heights. 


Figure  4,  Wind  perturbation  spectra  from  Gardner’s  gravity  wave  model. 


Fourier  inversion  using  random  phase  is  applied  to  the  spectra  to  generate  realizations  of  wind 
perturbation  profiles.  A  wind  perturbation  profile  is  generated  for  each  of  the  five  spectra.  A 
composite  profile  is  then  computed  by  shading  each  profile  spatially  with  a  Gaussian  filter  and 
then  summing  them  together,  where  Gaussian  filter  half-power  points  are  set  to  the  midpoint 
between  each  of  the  spectral  heights.  To  model  range-dependent  variability,  a  dominant 
horizontal  length  scale  is  defined,  and  Gaussian  weighting  functions  are  used  to  combine  the 
wind  perturbation  profiles.  Figure  5  gives  an  example  realization  of  a  wind  perturbation  field 
generated  from  the  gravity  wave  spectral  mode  using  a  horizontal  correlation  length  of  500  km 
See  Section  4.6.4  for  an  example  application  of  using  the  gravity  wave  model  to  compute 
propagation  variability. 
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Hor.  velocity  perturbations  (m/s) 


Figure  5.  Example  realization  of  range-dependent  wind  perturbation  field. 


4.4  Network  Performance. 

The  enhancements  to  InfraMAP  include  the  capability  to  evaluate  network  performance.  This 
capability  is  accessed  through  use  of  a  new  Source  Localization  form  (Figure  6),  and  it  enables 
the  localization  of  events  based  upon  measurement  data  and  modeling  results.  See  Section  4.7.4 
for  an  example  application  of  the  network  performance  capabilities. 
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Figure  6.  Source  Localization  form  that  provides  the  interface  to  the  network  performance 
capabilities. 
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The  localization  algorithm  is  a  Maximum  Likelihood  Estimator  (MLE)  implementation  (Kay, 
1993)  in  geophysical  coordinates  (Perl,  1981),  The  inputs  to  the  MLE  algorithm  are  a  set  of 
back  azimuth  and  arrival  time  measurements  and  associated  uncertainties  from  a  finite  number  of 
spatially  separated  sensors.  The  uncertainties  can  include  both  measurement  error  and 
uncertainty  associated  with  the  propagation,  as  discussed  in  later  sections.  The  state  of  interest  is 
positioned  in  spherical  coordinates  (lat.  Ion)  and  origin  time.  The  basic  algorithm  can  be  applied 
to  subsets  of  the  measurements,  but  note  that  the  state  vector  does  not  include  an  origin  time  if 
only  azimuthal  information  is  used. 

The  measurement  equation  relates  the  measurements  to  the  state.  For  azimuth,  it  is  determined 
directly  from  the  spherical  geometry.  The  arrival  time  is  related  to  the  state  by  defining  the 
average  sound  velocity,  or  signal  velocity ,  along  the  propagation  path,  where  the  signal  velocity 
is  defined  as  the  total  great  circle  range  divided  by  the  travel  time  (measured  arrival  time  minus 
origin  time),  Because  the  measurement  equation  is  nonlinear,  it  must  be  expanded  around  a 
reference  state  and  a  Jacobian  matrix  computed.  The  Jacobian  matrix  is  the  partial  derivative  of 
the  measurement  equation  with  respect  to  the  state  variables. 

The  Jacobian  matrix  and  measurement  inputs  are  combined  to  form  the  Fisher  information 
matrix,  from  which  an  unbiased  estimate  of  the  state  can  be  made.  Implicit  in  this  estimate  is  the 
assumption  that  the  uncertainties  can  be  represented  as  zero  mean  Gaussian  distributions, 
Depending  on  the  accuracy  of  the  initial  reference  state,  it  may  be  necessary  to  iterate  until 
convergence  is  reached.  A  covariance  matrix  that  specifies  the  uncertainty  in  the  state  estimate 
is  also  computed.  It  is  a  non-trivial  but  straightforward  analytic  exercise  to  map  the  covariance 
matrix  into  an  uncertainty  ellipse  over  a  specified  confidence  bound  [Farrell  and  LePage,  1996]. 
The  confidence  bound  is  usually  chosen  as  plus/minus  2  standard  deviations  or  95%. 


4.5  Reports  to  the  Research  Community. 

Under  this  contract,  BBN  participated  in  the  22nd  Annual  Seismic  Research  Symposium  |  Gibson 
and  Norris,  2000],  23rd  Seismic  Research  Review  [Norris  and  Gibson,  2001a],  and  24th  Seismic 
Research  Review  [Norris  and  Gibson,  2002a],  The  participation  included  submission  of  papers 
for  the  proceedings,  presentation  of  posters,  and  attendance  at  meeting  activities. 

BBN  also  has  given  technical  reports  in  several  other  meetings.  We  presented  research  on 
infrasonic  propagation  from  bolides  at  the  LANL  sponsored  Superbolide  Workshop  [Norris  et 
al. ,  2001].  Our  work  on  bolide  propagation,  including  an  analysis  of  elevated  ducts,  was  given  at 
the  Infrasound  Technology  Workshop  Kauila-Kona  [Norris  and  Gibson,  2001b],  At  the 
Infrasound  Technology  Workshop  in  De  Bilt,  we  presented  work  on  environmental  variability 
effects  [Gibson  and  Norris,  2002a].  Finally,  our  work  on  localization  accuracy  of  the  Pacific 
bolide  was  given  during  a  special  session  on  infrasound  at  the  Acoustical  Society  of  America 
meeting  in  Cancun  [Norris  and  Gibson,  2002b], 
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4.6  Sensitivity  Analyses. 

Several  sensitivity  studies  have  been  carried  out  under  this  contract.  They  have  focused  on 
applying  the  enhanced  environmental  and  propagation  capabilities  to  evaluate  the  effect  of 
various  environmental  variables  and  conditions  on  infrasonic  propagation. 

4.6.1  Environmental  Time  Evolution. 

A  new  feature  in  InfraMAP  is  the  ability  to  allow  for  time  evolution  of  the  environment  during  a 
propagation  run.  This  capability  can  be  used  with  both  ray  tracing  and  PE  codes.  Alternatively, 
a  user  can  freeze  the  time  evolution  at  a  fixed  time  (e.g.,  the  event  time).  With  time  evolution 
enabled,  the  climatological  profiles  at  a  point  along  the  propagation  path  are  determined  based 
on  the  origin  time  plus  the  elapsed  time  of  the  ray  (or  of  the  step  along  the  PE  slice).  In  a 
“frozen”  time  scenario,  all  climatological  profiles  are  referenced  to  the  same  time. 

A  study  was  conducted  on  the  significance  of  the  diurnal  effect  on  propagation  by  comparing  ray 
and  PE  results  both  with  and  without  time  evolution  enabled.  In  this  study,  the  source  event  time 
was  used  in  the  case  where  time  evolution  was  frozen,  The  HWM-93  and  MSISE-90 
climatologies  were  used  to  model  diurnal  effects.  Because  diurnal  effects  are  most  prevalent  in 
the  thermosphere,  thermospheric  paths  were  selected  for  analysis. 

The  scenario  chosen  for  modeling  is  as  follows.  The  source  was  located  in  central  Asia,  at  78° 

E,  50.5°  N.  The  source  time  was  15 -Sept  (Day  258)  at  10  UT.  A  propagation  path  of  2500  km 
to  the  north  was  modeled  using  both  eigenrays  and  a  PE  solution.  The  path  corresponds  to  a 
travel  time  of  approximately  three  hours  (10  -  13  UT).  Based  on  the  climatology,  both  the  zonal 
and  meridional  winds  shift  direction  in  the  thermosphere  during  this  time  window,  over  the 
region  of  infrasound  propagation.  The  nominal  propagation  path  is  shown  in  Figure  7. 
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Figure  7.  Nominal  propagation  path  (arrow)  for  time  evolution  sensitivity  study. 


The  effective  sound  speed  in  the  direction  of  propagation  along  the  great  circle  path  is  shown  in 
Figure  8.  The  panel  on  the  left  shows  the  evolving  environment,  and  the  panel  on  the  right 
shows  the  frozen  environment.  Based  on  the  climatology,  time  evolution  increases  the  turning 
heights  of  infrasound  in  the  thermospheric  region. 


Range  (km)  Range  (km) 

Figure  8.  Effective  sound  speed  along  great  circle  path  for  evolving  (left  panel)  and  frozen 

(right  panel)  environment. 


Parabolic  Equation.  The  results  predicted  by  the  PE  model  in  these  two  environmental 
scenarios  are  shown  in  Figure  9.  Again,  the  panel  on  the  left  shows  the  evolving  environment. 
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and  the  panel  on  the  right  shows  the  frozen  environment,  The  effect  of  the  environment  on 
turning  height  can  be  seen  by  examination  of  the  amplitude  fields  in  the  region  of  100  to  120  km 
altitude.  The  predicted  focusing  effects  at  the  ground  also  differ  somewhat  in  the  two  scenarios. 


Range  (km)  Range  (km) 


Figure  9.  PE  model  results  for  evol  ving  (left  panel)  and  frozen  (right  panel)  environment. 

InfraMAP's  capabilities  can  be  used  to  examine  the  effects  at  the  ground  in  more  detail.  In 
Figure  1 0,  amplitude,  relative  to  the  source,  is  shown  at  the  ground  (0  km  altitude)  for  both 
evolving  and  frozen  environments.  At  certain  ranges  (e,g.,  800  km,  1  800  km),  amplitude 
differences  on  the  order  of  20  dB  can  be  seen. 
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Figure  10.  PE  model  predictions  along  ground  for  evolving  and  frozen  environment. 
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Ray  Tracing.  When  ray  tracing  was  applied  to  the  two  scenarios,  three  thermospheric  eigenrays 
were  identified  for  each  case.  They  are  shown  in  Figure  1 1  for  both  the  evolving  and  frozen 
environments.  The  travel  time  differences  between  the  two  sets  of  rays  were  small,  on  the  order 
of  0.2%,  and  the  azimuth  deviation  differences  were  0.3  degree  or  less.  However,  there  are 
significant  differences  in  the  thermospheric  turning  heights  (up  to  approximately  5  km 
differences),  in  a  region  of  the  atmosphere  where  absorption  is  high  and  increasing  rapidly  with 
altitude. 


500  1000  1500  2000  2500 

Range  (km) 

Figure  1 1 .  Ray  tracing  paths  for  evolving  and  frozen  environment. 


Acoustic  absorption  is  calculated,  using  the  classical  absorption  model,  for  each  of  the  rays 
shown  above,  and  presented  in  Figure  12.  Large  differences  in  absorption  can  be  seen  between 
the  two  cases  for  those  paths  that  go  high  into  the  thermosphere  (the  tw:o  uppermost  curves).  The 
integrated  absorption  along  the  ray  path  is  important  because  it  determines  the  observability  of 
the  energy  at  a  receiver  location. 
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Ray  Trace:  Attenuation 
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Figure  12.  Attenuation  along  ray  path  for  evolving  and  frozen  environment. 


4.6.2  Atmospheric  Absorption. 

In  this  study,  the  integrated  absorption  mode!  [Sutherland  and  Bass,  1996]  is  applied  to 
characterize  the  attenuation  of  infrasound  at  various  frequencies.  Figure  1 3  gives  example 
profiles  of  predicted  absorption  coefficients  at  0.2  and  0,5  Hz.  Translation  and  diffusion  losses 
are  plotted  together  as  classical  loss,  of  which  the  diffusion  loss  makes  up  less  than  1%  of  the 
sum. 
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Sutherland  Attenuation  Model  Results  for  Frequency  =  Q.2Hz  Sutherland  Attenuation  Model  Results  for  Frequency  =  0.5Hz 


Figure  13.  Absorption  coefficients  computed  from  a  low  frequency,  high  altitude  model 


[Sutherland  and  Bass,  1996]  at  0.2  and  0.5  Hz. 

Below  60  km,  vibrational  loss  dominates,  but  the  total  absorption  is  still  relatively  small. 
Propagation  of  a  0.5  Hz  signal  over  a  500  km  path  would  only  experience  about  0.5  dB  of 
attenuation  in  this  region  (neglecting  spreading  loss).  The  effects  are  much  different  within  the 
thermosphere.  At  120  km,  the  absorption  coefficient  increases  dramatically  with  frequency;  at 
0,2  Hz,  it  is  0. 1  dB/km,  and  at  0.5  Hz  it  increases  to  1  dB/km.  Even  acoustic  signals  that  sample 
these  regions  for  a  small  percentage  of  their  propagation  path  can  experience  significant 
attenuation.  As  an  example,  a  0.5  Hz  signal  that  spends  10%  of  its  500  km  path  at  120  km  or 
above  would  have  an  absorption  loss  of  at  least  50  dB. 

Comparisons  were  made  of  PE  amplitudes  with  and  without  the  atmospheric  absorption  effects. 
Results  are  shown  in  Figure  14  for  a  source  at  the  ground.  The  amplitude  predictions  are  made 
for  upwind  conditions,  so  a  stratospheric  return,  with  an  upper  turning  height  of  40-60  km.  is  not 
present.  In  the  cases  where  absorption  is  not  included,  significant  thermospheric  energy  reaches 
a  receiver  on  the  ground  at  500  km  range,  With  absorption  included,  a  received  signal  at  500  km 
range  would  still  contain  a  thermospheric  return  at  0.2  Hz,  although  the  amplitude  is  reduced. 
However,  at  0.5  Hz,  the  signal  experiences  strong  absorption  at  altitudes  between  1 10  and  120 
km,  and  any  thermospheric  energy  arriving  at  the  receiver  is  heavily  attenuated. 
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Figure  14.  Comparison  of  PE  amplitude  predictions  with  and  without  absorption. 


Absorption  effects  at  0.5  Hz  are  apparent  in  the  thermosphere,  where  energy 


propagating  above  1 1 5  km  is  heavily  attenuated. 


Infrasonic  observations  are  ordinarily  made  on  the  ground.  Thus,  it  is  of  interest  to  look  at  the 
frequency  dependence  of  PE  amplitude  prediction  at  ground  level  (0  km  altitude).  Figure  1 5 
gives  PE  amplitude  at  the  ground  along  a  2600  km  path  and  a  broad  frequency  band.  The  figure 
suggests: 

•  A  trend  of  increasing  frequency  absorption  with  range 

•  Evolution  of  sharp  shadow  zone  boundaries  as  frequency  increases 

•  Shadow  zone  formation  over  limited  frequency  bands 

It  is  important  to  note  that  scattering  from  atmospheric  turbulence  will  tend  to  smear  energy  path 
boundaries.  As  a  result,  the  deep  attenuation  predicted  in  the  shadow  zones  would  be  reduced. 
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Figure  15.  PE  amplitude  predictions  along  the  ground  over  a  frequency  band  from  0.01  to 
2  Hz. 

4.6.3  Global  Duct  Heights. 

In  some  recent  studies  of  infrasonic  events  by  the  research  community,  disagreement  between 
predictions  and  observations  has  suggested  the  importance  of  elevated  duct  propagation.  There 
have  been  certain  observations  that  are  characterized  by  signal  velocities  typical  of  stratospheric 
ducting  that  cannot  be  modeled  with  stratospheric  ducts.  It  is  hypothesized  that  stratospheric 
energy  in  an  elevated  duct  is  observed  at  the  ground  through  diffraction  and  scattering.  This 
study  is  focused  on  quantifying  the  distribution  of  elevated  stratospheric  ducts. 

These  computations  are  based  on  using  HWM-93  and  MSIS-90  environmental  models  on 
January  1,  2001  at  12  UT.  The  date  corresponds  to  a  season  where  strong  zonal  winds  are 
present.  Stratospheric  ducting  would  be  less  pronounced  during  the  transitional  periods  in  the 
spring  and  fall  when  the  zonal  winds  change  direction. 

The  stratospheric  ducts  are  evaluated  along  the  maximum  wind  direction.  Sound  speed  profiles 
are  computed  for  the  no  wind  (static)  and  counter-wind  (effective)  cases.  The  counter-wind 
direction  is  chosen  to  evaluate  the  presence  of  elevated  stratospheric  ducts  for  the  case  where  no 
ground-based  stratospheric  ducts  exist. 

Elevated  stratospheric  ducts  have  the  same  upper  boundary  as  stratospheric  ducts  but  have  a 
lower  boundary  that  is  above  the  earth's  surface.  The  presence  of  elevated  stratospheric  ducts  is 
found  by  computing  the  height  of  ducts  starting  at  a  given  lower  boundary.  The  lower  boundary 
heights  evaluated  here  are  5  and  10  km.  The  presence  of  stratospheric  ducts  was  also  found 
using  a  lower  height  ot  0  km.  When  the  duct  height  is  above  1 00  km.  no  stratospheric  duct  was 
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found  and  only  a  thermospheric  duct  exists.  Figure  16  shows  the  duct  heights  for  the  10  km 
lower  boundary. 


Figure  16.  Duct  heights  for  lower  boundary  of  10  km  modeled  for  1  Jan  01. 


The  duct  height  features  in  Figure  16  can  be  divided  into  three  general  regions:  N,  C,  and  S.  The 
representative  profiles  for  each  region,  using  the  10  km  lower  boundary,  are  given  in  Figure  1 7 
to  19.  Two  sound  speed  profiles  are  shown:  the  sialic  case  with  no  wind,  and  the  counter- wind 
case.  The  green  arrow  shows  the  elevated  duct  that  forms  for  the  counter-wind  profiles.  In  the 
southern  region  (S),  no  elevated  stratospheric  ducts  exist  as  the  upper  height  boundary  is  in  the 
thermosphere.  For  the  central  region,  elevated  stratospheric  ducts  exist  with  upper  height 
boundaries  of  approximately  50  km.  In  the  northern  region,  elevated  stratospheric  ducts  form 
that  have  upper  height  boundaries  below  30  km. 


sound  speed  (m/s) 


Figure  17.  Static  and  counter- wind  Sound  speed  profiles  in  region  S. 
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sound  speed  (m/s) 

Figure  18.  Static  and  counter-wind  Sound  speed  profiles  in  region  C. 


sound  speed  (m/s) 

Figure  19.  Static  and  counter-wind  Sound  speed  profiles  in  region  N. 

Similar  calculations  were  done  at  lower  duct  heights  of  0  and  5  km.  At  5  km.  only  a  small  band 
of  elevated  stratospheric  ducts  exist.  This  band  is  centered  at  approximately  20  deg  North  in 
latitude. 

This  study  suggests  that  conventional  propagation  analysis  would  not  predict  any  stratospheric 
propagation  paths  in  the  counter  wind  direction.  However,  by  considering  elevated  stratospheric 
ducts  with  lower  heights  above  the  ground,  stratospheric  propagation  can  be  predicted  in  certain 
regions  of  the  globe,  with  the  size  of  the  region  dependent  on  lower  boundary  height.  To  use 
these  predictions  with  ground-based  measurements,  an  assumption  must  be  made  as  to  how 
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energy  leaks  out  of  the  duct  and  down  to  the  ground.  Diffraction  and  scattering  are  the  most 
likely  mechanisms. 

4.6.4  Propagation  Variability. 

The  goal  of  the  propagation  variability  study  is  to  quantify  the  bounds  in  travel  time  and  azimuth 
that  can  be  expected  as  driven  by  the  environmental  variability.  To  perform  the  study,  multiple 
wind  perturbation  realizations  are  generated,  and  a  Monte  Carlo  simulation  is  executed  where 
multiple  rays  are  traced  through  the  sum  of  mean  and  perturbed  environmental  fields.  The  two 
ray  parameters  (travel  time,  azimuthal  deviation)  are  calculated  for  each  perturbation.  The 
sensitivity  of  ray  tracing  calculations  to  variability  in  wind  profiles  is  then  quantified  by 
computing  the  mean  and  variance  of  the  predicted  distributions. 

Propagation  variability  is  studied  along  the  1700  km  path  from  the  23  April  2001  Pacific  bolide 
to  IS57  (Pinon  Flat,  California).  The  reference  ray  propagates  in  an  elevated  stratospheric  duct 
with  a  lower  boundary  height  of  approximately  8  km.  In  the  first  perturbation  study,  the 
environmental  variability  is  modeled  using  the  range-independent  gravity  wave  model.  Thus,  a 
signal  perturbation  profile  is  used  at  all  ranges  for  each  realization.  Twenty  perturbed  rays  are 
computed  and  shown  in  Figure  20  along  with  the  reference  ray.  The  perturbed  rays  deviate  from 
the  reference  ray  but  still  follow  the  same  general  ducted  path. 
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Figure  20.  Reference  and  perturbed  rays  using  the  range-independent  gravity  wave  model. 
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The  travel-time  distribution  for  the  perturbed  rays  is  given  in  Figure  21 .  The  small  sample  size 
tends  to  distort  the  distribution  but  it  roughly  resembles  a  normal  distribution.  The  mean  travel 
time  is  581 3  s  and  the  standard  deviation  is  44  s.  The  minimum  and  maximum  are 
approximately  100  s  plus  and  minus  the  mean  value. 
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Figure  21.  Travel-time  distribution  using  the  range-independent  gravity  wave  model. 


The  same  variability  is  also  studied  using  a  range-dependent  gravity  wave  model.  Here  the 
wind  perturbation  realizations  are  not  profiles  but  range-height  fields,  as  shown  in  Figure  5. 

Two  perturbed  rays  along  with  the  reference  ray  are  shown  in  Figure  22.  These  perturbed  rays 
deviate  fundamentally  from  the  reference  ray  in  that  they  are  trapped  in  secondary  transient 
ducts.  These  transient  ducts  result  from  the  wind  perturbation  field  and  exist  over  the  horizontal 
length  scale  of  the  gravity  wave  model  (selected  to  be  500  km  in  this  case).  Perturbed  rays  that 
are  trapped  in  transient  ducts  make  up  about  15  percent  of  the  total  number  of  realizations. 
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Figure  22.  Reference  and  two  perturbed  rays  using  the  range-dependent  gravity  wave  model. 


The  associated  travel  time  distribution  is  shown  in  Figure  23.  The  mean  is  5788  s  and  the 
standard  deviation  doubles  to  85  s.  Of  greatest  interest  is  that  the  distribution  now  has  outliers 
outside  of  the  main  distributions,  corresponding  to  the  perturbed  rays  trapped  in  the  transient 
ducts.  The  resulting  travel  times  are  approximately  300  s  earlier  than  the  mean,  and  the 
corresponding  signal  velocities  are  in  the  range  of  .3 10  km/s.  These  velocities  are  faster  than  the 
mean  value  of  0.294  km/s,  the  latter  more  typical  of  stratospheric  rays. 
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Figure  23.  Travel-time  distribution  using  the  range-dependent  gravity  wave  model. 

This  study  has  shown  that  new  propagation  features  can  be  seen  with  the  introduction  of  range- 
dependence  in  the  environmental  variability  models.  The  range-dependent  gravity  wave  model 
more  realistically  models  the  horizontal  scales  associated  with  gravity  wave  observations.  The 
perturbed  rays  suggest  that  energy  can  get  trapped  within  shallow  transient  ducts  that  form 
between  gravity  wave  layers.  One  modeling  effect  of  this  propagation  mode  is  a  reduction  in  the 
overall  travel  time  and  increase  in  signal  velocity. 


4.7  Validation  Studies. 

Several  validation  studies  have  been  performed  using  infrasonic  sources  of  opportunity. 
Launches  of  the  space  shuttle  provide  a  strong  moving  source  that  travels  from  the  ground  into 
the  thermosphere.  Bolides  are  another  excellent  infrasonic  source.  They  move  at  hypersonic 
speeds  and  fragment  (explode)  at  some  point  along  their  trajectory. 

4.7.1  Space  Shuttle. 

Rocket  launches  may  serve  as  useful  ground  truth  data  for  infrasound  [McLaughlin  el  al.,  2000] 
and  also  represent  an  excellent  source  of  opportunity  for  model  validation.  Space  shuttle 
launches  have  been  observed  at  Los  Alamos  (DL1AR)  and  at  Lac  du  Bonnet  (IS  10),  among  other 
locations. 
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Rocket  and  missile  launches  from  the  Eastern  US  were  observed  extensively  in  the  1960’s  and 
1970's  at  Palisades,  NY  and  elsewhere,  and  a  series  of  reports  describing  them  (for  example, 
[Balachandran  and  Donn,  1971])  were  issued  by  scientists  at  Lamont-Doherty  Geological 
Observatory  and  the  US  Army  Electronics  Command.  A  number  of  important  findings  resulted 
from  this  early  work,  including: 

•  Characterization  of  seasonal  dependency; 

•  Quantification  of  wind  effects  on  azimuth; 

•  Identification  of  two  distinct  source  regions,  one  near  the  launch  site  and  one  near  the  re¬ 
entry  location  of  the  first  stage; 

•  Association  of  source  location  with  high-velocity  portions  of  the  trajectory  in  high- 
density  regions  of  the  atmosphere. 

NASA’s  space  shuttle  is  launched  from  Cape  Kennedy,  Florida  a  few  times  per  year.  During  the 
spring  and  summer  of  1999,  the  launches  of  two  shuttle  missions  were  observed  at  infrasound 
stations  DLIAR  and/or  IS1 0. 

BBN  used  a  commercially  available  launch  simulator  software  package  to  determine  the 
trajectories  for  these  specific  missions,  using  actual  launch  parameters.  The  software  was  also 
used  to  estimate  trajectories  for  the  solid  rocket  boosters,  which  are  released  from  the  orbiter 
approximately  2  minutes  into  the  ascent.  Modeled  trajectories  for  the  two  observed  shuttle 
missions,  STS-96  (27  May  99)  and  STS-93  (23  July  99)  are  shown  in  Figure  24.  Launch  ascent 
trajectories  for  the  STS-96  orbiter  and  solid  rocket  boosters  (SRBs)  are  shown  in  Figure  25.  (The 
modeled  STS-93  ascent  trajectories  are  very  similar  to  those  for  STS-96.) 
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Figure  24.  Space  shuttle  trajectories  for  two  missions. 
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STS-96  Launch  Ascent  Trajectory  (27  May  99) 


InfraMAP  was  used  to  determine  eigenrays  to  DLIAR  from  points  along  the  STS-96  trajectory, 
with  HWM-93  and  MSISE-90  used  for  environmental  characterization.  The  continuously 
moving  source  was  modeled  as  a  series  of  discrete  sources  separated  in  space  and  time.  A  source 
was  modeled  every  10  seconds  from  the  launch  time  out  to  5  minutes  after  launch  for  the  orbiter 
and  from  200  seconds  out  to  6  minutes  after  launch  for  the  solid  rocket  boosters  (SRB).  For  each 
eigenray,  an  arrival  azimuth  and  an  arrival  time  (referenced  to  the  launch  time)  were  determined. 
Results  are  shown  in  Figure  26.  Stratospheric  rays  and  thermospheric  rays  are  depicted 
separately  for  both  the  orbiter  and  the  SRB.  Also  shown  in  the  figure  (as  red  asterisks)  are 
results  from  the  observation  at  DLIAR.  determined  by  analysis  using  InfraTool  and  MatSeis 
software,  and  provided  by  Dr.  Rod  Whitaker  of  DoE/LANL  and  Dr.  David  Brown  of  CMR. 
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STS- 96  to  DLIAR 


time  r,e.  launch  (sec) 


Figure  26.  Predicted  and  observed  infrasound  arrivals  from  STS-96  to  DLIAR. 

The  two  primary  observed  arrivals  (at  approximately  103.5  and  99  degrees)  are  reasonably  well 
modeled  by  the  stratospheric  rays  from  the  orbiter  (blue  circles)  and  boosters  (green  squares), 
respectively.  However,  a  bias  in  azimuth  of  approximately  2  degrees  can  be  seen  for  both 
primary  arrivals.  Further  modeling  of  the  launch  event  using  an  updated  atmospheric 
characterization,  in  order  to  see  if  travel  time  and  azimuth  predictions  could  be  improved,  would 
be  of  great  interest  in  this  case. 

Infrasound  from  space  shuttle  launches  has  also  been  observed  at  an  infrasound  array  operated 
by  the  Army  Research  Lab  (ARI.)  at  Blossom  Point,  Maryland.  [Tenney  and  Noble,  2000 j 
InfraMAP  was  used  to  model  propagation  from  the  launch  of  the  STS-88  mission  (04  Dec  98), 
again  using  HWM-93  and  MSISE-90.  The  STS-88  trajectory  was  similar  to  that  of  STS-96 
shown  above.  Similar  time  steps  and  durations  were  used  for  the  source  model.  A  spectrogram 
of  the  observation  at  Blossom  point  (data  provided  by  Mr.  Steve  Tenney  and  Dr.  John  Noble  of 
ARL)  is  shown  in  Figure  27.  An  azimuth  vs.  time  diagram  is  shown  in  Figure  28  The  white 
“dotted  line”  boxes  have  been  added  to  Figure  28  to  denote  the  primary  arrivals. 
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Figure  27.  Observed  infrasound  signal  from  STS-88  to  Blossom  Point. 
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Figure  28.  Observed  infrasound  arrivals  from  STS-88  to  Blossom  Point. 


39 


Eigenray  results  from  InfraMAP  are  shown  in  Figure  29.  Overlaid  on  the  predictions  are  red 
“dotted  line”  boxes  that  correspond  to  the  arrivals  as  depicted  in  Figure  28.  In  this  case,  the 
increasing  azimuth  trend  is  well  modeled  by  both  the  stratospheric  and  thermospheric  rays. 

There  is  a  bias  in  travel  time  of  approximately  300  seconds  between  the  observed  arrivals  and  the 
modeled  stratospheric  rays.  The  modeled  stratospheric  rays  from  the  boosters  (green  squares) 
arrive  before  those  from  the  orbiter  (blue  circles)  because  the  source  is  moving  supersonically 
toward  the  array. 


STS- 88 


Figure  29.  Predicted  and  observed  infrasound  arrivals  from  STS-88  to  Blossom  Point  (red 

boxes  depict  main  observed  arrivals,  as  shown  in  Figure  27). 


Further  research  into  several  technical  issues  has  been  ongoing  under  a  separate  contract  effort, 
including:  modeling  improvements  achievable  with  near- real -time  updates  as  compared  to 
HWM-93  and  MSISE-90;  further  understanding  of  the  infrasound  source  mechanism  in  order  to 
identify  the  regions  of  the  trajectory  (altitude,  velocity,  etc.)  that  contribute  most  strongly  to 
observed  signals;  and  quantifying  attenuation  along  the  ray  paths  to  support  identification  of 
observed  phases. 
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4.7.2  El  Paso  Bolide 


On  9  Oct  1997,  a  large  bolide  traveled  above  Texas  near  El  Paso.  The  height  of  the  bolide,  as 
determined  from  satellite  observations,  was  29  km.  This  bolide  was  detected  at  both  the  Los 
Alamos  DLIAR  array  (range  445  km)  and  the  Southern  Methodist  University  TXIAR  array 
(range  359  km).  The  received  data  at  TXIAR  has  been  studied  to  evaluate  the  ability  to  verify  or 
refine  the  source  height  estimate. 

To  model  the  bolide  signal,  eigenrays  were  computed  over  a  range  of  source  height  from  0  to  50 
km.  Figure  30  shows  the  eigenray  solutions  at  30  km  source  height.  Over  the  range  of  heights, 
there  are  1  -3  stratospheric  eigenrays  and  2-3  thermospheric  eigenrays. 


Figure  30.  Eigenray  solution  between  El  Paso  bolide  at  source  height  of  30  km  and  receiver  at 


TXIAR. 

Synthetic  waveforms  were  next  computed  for  each  source  height  and  mapped  into  an  image  as  a 
function  of  arrival  time  and  source  height  (bottom  of  Figure  31).  A  Gaussian  waveform  of 
standard  deviation  10  sec  was  used  as  the  source  function.  For  sources  near  the  ground,  single 
eigenrays  can  be  observed,  which  then  split  into  multiple  eigenrays  at  the  height  increases.  The 
split  rays  for  a  given  v-shaped  branch  correspond  to  two  similar  paths,  one  propagating  up  at  the 
source,  the  other  propagating  down,  reflecting  off  the  ground,  and  then  propagating  up.  The 
eigenray  image  can  be  compared  to  the  measured  waveform  (top  of  Figure  31),  and  a  source 
height  estimate  can  be  made  by  correlating  the  measured  arrival  times  with  the  predicted 
amplitudes.  The  three  measured  arrivals  match  the  predictions  at  source  heights  of  33,  26,  and 
23  km,  respectively.  Although  there  is  a  spread  in  these  heights,  they  are  in  general  agreement 
with  the  satellite-determined  source  height  of  29  km. 
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Figure  31.  Comparison  of  measured  arrival  times  (above)  with  synthetic  waveform  prediction 


from  eigenrays  solved  over  0  to  50  km  (below). 


4.7.3  Tagish  Lake  Bolide 

On  1 8  January  2000,  sensors  aboard  DOD  satellites  detected  the  impact  of  a  meteoroid  at 
16:43:43  UT  near  Whitehorse  in  the  Yukon  Territory,  Canada.  The  object  detonated  at  an 
altitude  of  25  km  at  60,25  decrees  North  latitude,  134.65  degrees  West  longitude  [Brown  et  al. . 
2002] 

An  inffasonic  signal  was  detected  at  the  IMS  station  in  Lac  du  Bonnet,  Canada,  a  range  of  2641 
km  from  the  bolide.  A  received  waveform  is  shown  in  Figure  32.  This  waveform  is  filtered  over 
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the  frequency  band  of  1 .0  - 1.5  Hz.  The  main  arrival  appears  at  an  arrival  time  of  approximately 
8200  sec.  There  may  also  exist  a  secondary  arrival  at  1 1,300  sec. 
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Figure  32.  Filtered  waveform  measured  at  Lau  de  Bonnet  IMS  station  (blue),  and  synthetic 
waveform  generated  from  eigenrays  (red).  Amplitude  units  are  arbitrary. 

Measurements  across  the  3  elements  of  the  array  were  beam  formed  using  an  F-K  processing 
algorithm.  The  results  were  inconsistent  across  the  observed  waveform,  possibly  due  to  low 
signal-to-noise  ratios  or  loss  in  coherence  of  the  signals  across  the  array,  resulting  from 
turbulence.  The  strongest  coherent  cluster  in  the  F-K  slowness  plane  was  found  using  a  time 
window  from  8400  -  8500  sec  and  a  frequency  band  of  0.2  -  1.0  Hz.  The  back  azimuth  for  this 
cluster  was  306.5  deg. 

Eigenrays  were  computed  for  this  path  and  are  shown  in  Figure  33.  Both  stratospheric  and 
thermospheric  rays  are  present.  The  stratospheric  rays  have  back  azimuths  at  approximately  310 
deg,  and  the  thermospheric  at  approximately  31 1  deg.  These  predictions  differ  by  several 
degrees  from  the  measured  value  of  306.5  deg. 
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Figure  33.  Eigenrays  from  Tagish  Lake  bolide  to  Lac  du  Bonnet  IMS  station. 

Synthetic  waveforms  were  generated  from  the  eigenrays  for  comparison  with  the  measurements. 
The  source  function  was  a  Gaussian  waveform  with  a  standard  deviation  of  10  sec.  The 
predictions  are  overlaid  on  the  waveforms  in  Figure  32.  A  bias  of  approximately  500  sec  is 
apparent  in  the  arrival  times.  This  difference  may  be  due  to  discrepancy  between  the  predictions 
of  temperature  and  wind  from  the  HWM-93  and  MSISE-90  models  and  the  actual  environmental 
conditions  during  the  propagation.  It  may  also  result  from  inadequate  modeling  of  the  source 
function. 

4.7.4  Pacific  Bolide 

In  this  study,  the  23  April  2001  Pacific  bolide  is  used  to  evaluate  the  network  localization 
capabilities  integrated  into  InfraMAP,  This  bolide  is  an  excellent  source  of  opportunity  because 
ground  truth  in  latitude,  longitude,  and  height  is  available  and  infrasound  from  the  bolide  was 
observed  at  a  large  number  of  stations. 

Figure  34  shows  the  location  of  the  bolide  and  the  six  infrasound  stations  with  observations  used 
to  localize  it.  Note  that  the  IS26  station  in  Freyung,  Germany  is  not  shown  in  the  figure,  the 
great  circle  path  length  for  this  case  being  over  9000  km 
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Figure  34.  Location  of  the  23  April  2003  Pacific  bolide  along  with  the  six  stations  used  to 
localize  it. 


The  localization  is  based  on  measured  arrival  times,  back  azimuths,  and  associated  measurement 
errors  recorded  at  the  stations.  A  baseline  localization  has  been  computed  using  these  data  and 
signal  velocity  tables  [Brown  et  al,  2001].  To  fully  utilize  the  integrated  localization 
capabilities  in  InfraMAP,  additional  estimates  of  azimuthal  deviation,  and  propagation 
uncertainty  due  to  the  environment  are  needed.  These  values  are  computed  by  using  an  event- 
specific  propagation  modeling  to  classify  the  ray  phase  (propagation  mode)  from  the  bolide  to 
each  station  [Norris,  2003;  Norris  and  Gibson,  2002a],  This  approach  also  estimates  signal 
velocity,  eliminating  the  need  for  look-up  tables. 

The  model-based  approach  starts  by  using  the  baseline  localization  as  an  initial  source  reference 
point  for  the  modeling.  A  fan  of  rays  is  reverse  propagated  from  each  station  back  toward  the 
source.  These  rays  are  then  identified  and  grouped  into  phases  (either  stratospheric  or 
thermospheric).  Two  ray  parameters  are  computed  for  each  phase:  signal  velocity  and  azimuthal 
deviation.  Likewise,  these  parameters  are  computed  using  the  reference  source  position  and 
station  measurements.  The  two  versions  of  the  parameters  arc  compared,  and  the  modeled  phase 
that  provides  the  best  match  is  identified. 

Total  azimuth  and  travel-time  uncertainty  is  computed  as  the  sum  of  the  measurement  error  plus 
a  contribution  based  on  the  uncertainty  in  characterizing  the  propagation.  The  latter  is  found  by- 
propagating  rays  through  perturbed  realizations  of  the  environment  and  building  up  a  set  of 
azimuth  and  travel-time  predictions  for  which  variances  can  be  extracted.  The  environmental 
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realizations  are  found  using  the  spectral  gravity  wave  model  that  models  the  relevant  spatial 
variability  in  winds. 

The  model -based  localization  is  calculated  using  the  same  MLE  algorithm  as  used  in  the  baseline 
localization  but  with  several  key  enhancements: 

•  the  mean  signal  velocity  for  the  modeled  phase  is  used; 

•  the  azimuthal  deviations  for  the  modeled  phase  is  used  and  added  to  the  measured  back 
azimuths  to  correct  for  horizontal  propagation  effects; 

•  new  estimates  of  travel  time  and  azimuthal  uncertainty  are  used  that  include  propagation 
variability  and  measurement  error. 

The  resulting  localizations  are  given  in  Figure  35.  A  localization  based  solely  on  azimuthal  data 
is  shown  in  yellow.  The  lack  of  complete  azimuthal  coverage  of  the  stations  is  evident  in  the 
elongation  of  the  error  ellipse  along  an  azimuth  aligned  with  most  of  the  stations.  A  localization 
using  only  arrival  time  data  is  shown  in  green.  Here  the  geographic  distributions  of  the  stations 
results  in  an  error  ellipse  aligned  perpendicular  to  the  general  station  direction.  Finally,  the 
localization  based  on  all  data  is  shown  in  cyan.  The  resulting  error  ellipse  is  significantly 
reduced  and  has  an  AOU  of  8700  square  km,  and  the  localization  has  a  small  29  km  miss 
distance  when  compared  to  ground  truth. 

Azimuthal  Equidistant  projection  centered  on  GT 


range  (km) 


Figure  35.  Model-base  localization  using  only  azimuth  data  (yellow),  travel-time  data  (green) 
and  all  data  (cyan). 
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The  model-based  approach  to  localizing  the  Pacific  bolide  showed  an  approximately  90% 
improvement  in  both  localization  accuracy  and  reduced  AOU  when  compared  with  the  baseline 
localization.  The  Pacific  bolide  example  illustrates  the  excellent  performance  of  the  model- 
based  localization  approach. 
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Section  5 

Conclusions  and  Recommendations 


Under  this  effort,  several  enhancements  have  been  made  to  the  InfraMAP  tool  kit  that  facilitate 
advanced  modeling  and  analysis  of  long-range  infrasonic  propagation.  These  enhancements 
have  been  applied  to  several  sensitivity  and  special  event  studies.  The  major  highlights  of  the 
project  are  detailed  below. 

Absorption 

An  absorption  model  developed  for  low  frequencies  and  high  altitudes  has  been  integrated  into 
the  package.  Absorption  is  particularly  significant  within  the  thermosphere.  PE  modeling  for 
upwind  propagation  over  a  500  km  path  indicates  that  thermospheric  energy  becomes  heavily 
attenuated  as  the  frequency  increases  from  0.2  to  0.5  Hz. 

Synthetic  Waveforms 

A  synthetic  waveform  capability  has  been  added  to  the  ray  model.  The  source  waveform  can  be 
specified  as  a  blast  waveform,  Gaussian  waveform,  or  user-defined  waveform.  The  source 
waveform  is  convolved  with  ray  arrival  times  and  weighted  by  ray  attenuation  factors  to 
synthesize  a  received  waveform. 

Environmental  Variability 

Gravity  waves  are  believed  to  be  the  dominant  source  of  environmental  variability  that  affects 
infrasonic  propagation.  A  gravity  wave  spectral  model  has  been  integrated  into  InfraMAP  based 
upon  scale-independent  diffusive  filtering  arguments.  From  the  model,  realizations  of  wind 
perturbation  profiles  are  generated.  The  profiles  are  based  upon  five  separate  spectra  in  the 
vertical  direction  and  one  dominant  Gaussian  length  scale  in  the  horizontal  direction. 

Network  performance 

A  network  performance  model  has  been  integrated  that  can  provide  source  localization  based  on 
measured  back  azimuth  and  arrival  times  across  a  local  network  of  infrasonic  stations.  The 
model  allows  for  integration  of  results  from  propagation  modeling,  including  the  effects  of 
environmental  variability. 

Propagation  variability 

Propagation  variability  has  been  studied  by  propagating  rays  through  perturbed  profiles,  and 
Monte  Carlo  statistics  found  for  travel  time  and  azimuthal  deviation.  These  characterizations 
have  been  used  to  evaluate  network  performance  over  a  group  of  sensors  by  computing  AOU 
confidence  bounds  for  the  source  localizations. 

Recent  variability  modeling  studies  indicate  the  presence  of  transient  ducts  that  can  significantly 
perturb  travel  times  and  signal  velocities.  These  ducts  form  around  localized  minima  in  the 
sound  speed,  and  their  extent  is  determined  by  the  horizontal  correlation  length  of  the  gravity 
waves. 
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El  Paso  Bolide 


Large  bolides  provide  a  strong  source  for  infrasound  and,  in  conjunction  with  ground  truth  data 
from  satellites,  have  been  used  to  evaluate  the  performance  of  the  models.  Source  height 
estimates  of  the  Ei  Paso  bolide  have  been  made  over  a  400  km  path.  Synthetic  waveforms  were 
computed  from  eigenray  solutions  over  a  range  of  hypothesized  heights  and  compared  with  the 
detected  arrivals  in  the  measured  waveform.  The  predicted  heights  all  fell  within  6  km  of  the 
ground  truth  height,  suggesting  that  multipath  arrivals  can  be  used  to  estimate  source  height. 

Pacific  bolide 

Data  from  the  Pacific  bolide  detected  at  six  stations  was  used  to  perform  multi-station  source 
localization.  The  localization  included  detailed  propagation  path  modeling  and  evaluation  of 
environmental  variability  effects.  The  resulting  location  provided  a  more  accurate  localization 
and  90  %  smaller  error  ellipse  w'hen  compared  to  a  table-based  localization  approach.  The  most 
important  factor  in  localization  accuracy  was  found  to  be  correct  classification  of  the  ray  paths. 

Elevated  duct 

Elevated  stratospheric  ducts  above  10  km  can  form  over  large  regions  of  the  globe.  These  ducts 
can  be  excited  by  either  elevated  sources  or  ground  source  energy  that  leaks  into  duct.  It  is 
hypothesized  that  this  energy  can  reach  ground-based  stations  through  scattering  and  diffraction. 

Recommendations  and  research  topics  that  would  support  improved  modeling  of  infrasonic 
events  are  listed  below. 

•  The  Tagish  Lake  bolide  predictions  suggest  the  possibility  that  a  bias  exists  in  the 
environmental  models  for  the  day  and  time  of  this  event.  Further  improvements  to  the 
environmental  models  in  InfraMAP  should  be  investigated.  In  addition,  data  assimilation 
of  empirical  models  with  in  situ  measurements  should  be  pursued  to  increase  the 
accuracy  of  the  propagation  predictions. 

•  The  gravity  wave  model  can  be  improved  by  including  the  geospatial  dependence  on  the 
amplitude  scale  factor.  Other  issues  that  should  be  examined  include:  mean  plus 
perturbation  superposition,  mean  flow  bias,  max  flow  limits,  anisotropy,  and  zonal  vs. 
meridional  flow  dependence. 

•  The  largest  errors  contributing  to  source  localization  uncertainty  are  generally  in  the  back 
azimuth  estimates.  Improvements  in  either  the  array  geometries  or  beam-forming 
algorithms  would  be  desirable  to  help  reduce  this  error  term. 

•  Network  performance  studies  need  to  be  done  over  a  variety  of  scenarios  and 
environmental  conditions.  In  areas  of  high  station  density,  the  localization  AOU  can  be 
small,  while  in  areas  of  sparse  coverage  they  can  become  quite  large.  Defining  these 
regions  of  best  and  w'orst  localization  accuracy  is  critical  to  evaluating  the  overall 
performance  of  an  infrasound  network, 
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•  General  analysis  of  future  bolides  needs  to  continue.  This  will  be  facilitated  by 
collaboration  with  the  satellite  community  and  with  the  international  community.  Bolide 
source  modeling  needs  to  improve,  particularly  in  the  areas  of  source  mechanics, 
transition  region  where  the  blast  wave  can  be  modeled  as  a  linear  propagating  wave,  and 
light  curve  to  acoustic  source  waveform  correlation. 

•  Ray  classification  techniques  need  to  be  evaluated  and  improved.  Use  of  additional 
discriminators,  such  as  phase  arrival  envelope  and  amplitude,  may  provide  new 
classification  approaches. 

•  Non-linear  propagation  effects,  such  as  waveform  steepening  and  multi-shock  merging, 
need  to  be  researched  for  infrasonic  sources.  These  effects  may  be  particularly 
significant  for  hypersonic  sources,  such  as  bolides  and  reentry  vehicles. 

•  The  frequency  dependence  on  wave  propagation  and  travel -time  predictions  needs  further 
research.  In  particular,  the  limits  into  the  applicability  of  ray  theory  to  infrasonic 
propagation  have  not  been  adequately  defined.  More  research  needs  to  be  done  in 
developing  and  evaluating  broadband  propagation  models,  such  as  the  Time- Domain 
Parabolic  Equation  (TDPE)  model. 
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Appendix 

Abbreviations  and  Acronyms 

AFTAC 

AGU 

CTBT 

GUI 

HARPA 

Air  Force  Technical  Applications  Center 

American  Geophysical  Union 

Comprehensive  Nuclear  Test  Ban  Treaty 

Graphical  User  Interface 

Hamiltonian  Ray-Tracing  Program  for  Acoustic  Waves  in  the 
Atmosphere 

HWM 

IMS 

InfraMAP 

LANL 

MSIS  or  MSISE 
NASA 

NOAA 

NPE 

NRL 

PE 

TDPE 

TXIAR 

UT 

WKJB 

Horizontal  Wind  Model 

International  Monitoring  System 

Infrasonic  Modeling  of  Atmospheric  Propagation 

Los  Alamos  National  Laboratory 

Extended  Mass  Spectrometer  -  Incoherent  Scatter  Radar  model 
National  Aeronautics  and  Space  Administration 

National  Oceanographic  and  Atmospheric  Administration 
Nonlinear  Progression  Wave  Equation 

Nava!  Research  Laboratory 

Parabolic  Equation 

Time  Domain  Parabolic  Equation 

Texas  Infrasound  Array 

Universal  Time 

Wentzel-Kramer-Brillouin  method 
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